arXiv: 1506.00526vl [astro-ph.HE] 1 Jun 2015 


Mon. Not. R. Astron. Soc. 000 , [[]-?? (2012) Printed 2 June 2015 


(MN IATeX style file v2.2) 


A viscous-convective instability in laminar Keplerian thin discs. II. 
Anelastic approximation. 


N. Shakura 1 *, K. Postnov 2,1 

1 Sternberg Astronomical Institute, Moscow M. V. Lomonosov State University, Universitetskij pr., 13, Moscow 119992, Russia 

2 Faculty of Physics, M. V. Lomonosov Moscow State University, Leninskie Gory, Moscow 119991, Russia 


Received ... Accepted ... 


ABSTRACT 

Using the anelastic approximation of linearised hydrodynamic equations, we investigate the 
development of axially symmetric small perturbations in thin Keplerian discs. The sixth-order 
dispersion equation is derived and numerically solved for different values of relevant physi¬ 
cal parameters (viscosity, heat conductivity, disc semi-thickness and vertical structure). The 
analysis reveals the appearance of two overstable modes which split out from the classical 
Rayleigh inertial modes in a wide range of the parameters in both ionized and neutral gases. 
These modes have a viscous-convective nature and can serve as a seed for turbulence in astro- 
physical discs even in the absence of magnetic fields. 
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1 INTRODUCTION 


In an attempt to understand hydrodynamic instabilities which can potentially initiate turbulence in accretion discs, in lShakura & Postnovl 
(.2015:) (Paper I) we have performed a local WKB-analysis of axially symmetric perturbations in thin accretion discs. It was found that 
under a special choice of wave vector direction of the perturbations, almost (but not completely) aligned with the disc symmetry plane, the 
presence of a microphysical viscosity, parametrized in terms of the mean-free path length to the disc scale ratio and taken into account in the 
dissipation function in the right-hand side of the energy equation, leads to the appearance of an overstable oscillating behaviour of one of two 
classical Rayleigh inertial modes. The instability was found in a wide range of perturbation wavelengths (expressed through the dimensional 
wavelength vector kr, where r is the disc radial scale), around kr ~ 30 - 100, in both fully ionized gases and neutral gases. The microphysical 
heat conductivity was taken into account through the dimensionless Prandtl number, Pr, which ranges from 0.052 for fully ionized plasma 
to 2/3 for neutral gases. We have found that the instability increment, reaching ~ 0.1 local Keplerian values, diminishes with decreasing the 
Prandtl number (e.g. due to the presence of a photon heat conductivity) and with increasing background vertical entropy gradient (expressed 
in terms of the Brunt-Vaisala frequency). Such a behaviour of the instability is in agreement with physically intuitive dumping effect of the 
heat conductivity and entropy gradients on the development of small radial perturbations propagating under a small angle to the disc plane. 

To make the physics as simple as possible, in Paper I we have used the Boussinesq approximation of the hydrodynamic equations, 
which assumes the incompressibility of the fluid in the continuity equation and neglects the Euler pressure variations in the energy equation. 
We argued that the incompressibility approximation is justified for radial perturbations with kr » 1, which may suggest that the discovered 
viscous instability of the inertial Rayleigh modes is real and not the result of the approximations used. 

In this paper we continue studying the viscous-convective instability in thin shear laminar flows found in our paper lshakura & Postnovl 
d2015h . Here we treat the problem in the anelastic approximation and take into account the vertical boundary conditions in the thin Keplerian 
discs. The anelastic approximation is the next approximation to the full system of hydrodynamic equations, but in which the term dp/dt = 0, 
i.e. the continuity equation takes the form div(pu ) = 0, which allows one to filter out sound waves. 

When considering sound-proof stratified flows, the use of the anelastic approximation is known to have some subtleties (see, for example, 
the recent analysis by I Vasil et al.1120131) and references therein). Special attention should be given to the energy equation, since the standard 
anelastic set of equations operates with adiabatic perturbations lOgura & Phillipsl d 19621) . Our analysis, in contrast, is heavily based on the 
viscous energy generation in the sheared flows, therefore the rigorous proof of the applicability of the anelastic approximation in this case is 
to be found. As a justification of this treatment we heuristically use the criterion that the linearised equations should not give rise to spurious 
modes with unphysical behaviour (e.g., unstable modes in the steady-state solid-body rotation case). 

Although our analysis is applicable for any sheared axially symmetric flow, we will be mostly concerned with thin Keplerian discs 
which have a wide range of phenomenological applications. This means that in the continuity equation of importance becomes the term 
~ l/p 0 (dpo/dz), which can be quite significant in thin discs and which we have neglected in the Boussinesq approximation. The vertical 
boundary conditions in thin discs are taken into account by solving the Sturm-Liouville problem for the "-part of perturbations. 

The main result of the paper is the dispersion equation Eq. {48}, which is a sixth-order algebraic equation for small perturbations in 
the form /(z)exp(iwf - k r r)). The solution of this equation signals the appearance of an overstable behaviour for two modes with the same 
negative imaginary part and real parts with equal absolute values but different signs which split from two classical inertial Rayleigh modes, 
in a wide range of k r r » 1. 

The structure of the paper is as follows. In Section [2] we write down the basic equations. In Section [3] we linearise the full system of 
equation in the anelastic approximation. We proceed with the derivation of the dispersion equation in Section[4] following by its numerical 
analysis in Section[5] Section [6] summarizes our findings. 


2 BASIC EQUATIONS 

The system of hydrodynamic equations reads: 


(i) mass conservation equation 


dp_ 

dt 


+ V • (pu ) = 0, 


0 ) 


In cylindrical coordinates for axially symmetric flows: 


V-(p«) 


1 d(pru r ) d(pu z ) 
r dr + dz 


( 2 ) 


(ii) Navier-Stokes equation including gravity force 

^ + (mV) • u = --Vp - VA, + N. (3) 

dt p 

Here i p g = -GM/r is the Newtonian gravitational potential of the central body with mass M , N is the viscous force. In cylindrical coordinates 
for axially symmetric flows: 


du r du r du r 

- + U r - + U 7 - 

dt dr ' dz 


Jh-L d JL + K , 

dr p dr 


(4) 



duj, du$ UfU^ 

dr z dz r 




( 5 ) 
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du- 


du. 


du. 


1 dp 




dt 


dr 


dz 


dz p dz 


The linearised viscous force components are specified in Appendix A of Paper I. 
(iii) energy equation 


( 6 ) 


pKT 

P 


ds 

~di 


+ (mV) • s 


= <2visc - V • F . 


(V) 


where 5 is the specific entropy per particle, <2visc is the viscous dissipation rate per unit volume, H is the universal gas constant, p is the 
molecular weight, T is the temperature, and terms on the right stand for the viscous energy production and the heat conductivity energy flux 
F, respectively. The energy flux due to the heat conductivity is 


V • F = V(-kVT) = -kAT -Wk-VT. ( 8 ) 

Note that both electrons and photons, and at low temperatures neutral atoms, can contribute to the heat conductivity (see Section[5]below). 
(iv) equation of state 

The equation of state for a perfect gas is convenient to write in the form: 

p = Ke s,cv p\ (9) 


where A' is a constant, Cy = l/(y - 1) is the specific volume heat capacity and y = c p tcv is the adiabatic index (5/3 for the monoatomic gas). 
We will also use the equation of state in the form 



( 10 ) 


where p is the molecular weight. 


3 LINEARISED EQUATIONS IN ANELASTIC APPROXIMATION 

The perturbed hydrodynamic variables can be written in the form x = Xq + Xi , where xo stand for the unperturbed quantities and 
x \ = (pi, Pi, u r% \, u z j , ) are small perturbations. In contrast to Paper I in which we considered the local WKB approximation, i.e. small per¬ 

turbations of density, pressure and velocity in the form x\(t, z, r) oc exp(\col-\k r r-ik z z), here we will take them in the form oc f{z) exp(iajt-ik r r) 
with the boundary conditions f(zo) = 0, f(-Zo) = 0, where zo is the disc semi-thickness. We will consider thin discs with zo/>' ~ uju ^0 <K 1 
(u s is the sound velocity). Below we shall omit subscript 1 for small perturbations of the velocity, unless stated otherwise. 

In this Section we will formultae the so-called anelastic approximation of h ydrodyna mic eq uations in which the sound wave perturba¬ 
tions are neglected by omitting the term dp/dt in the continuity equation lOeura & Phillips! 11962h . 

The linearised hydrodynamic equations are written as follows. 


(i) Continuity equation 

The anelastic approximation for gas velocity u is V • p 0 u = 0: 


du z 


1 r)n,. 




(ii) Dynamic equations 

The radial, azimuthal and vertical components of the Navier-Stokes momentum equation are, respectively: 


i cju r - 20.u,/, = i k, 


Pi_ Pi dp 0 
Po P 5 dr 


vkz 


u r + v- 


d 2 u r 

dz 2 


(ID 


( 12 ) 


K ,0 

, H- U, — —VkrUp, 

20 r * 


+ V' 


d 2 u t f, 

dz 2 


(13) 


1 dp 1 Pi dp 0 

p dz + pi dz 


,2 ^ d 2 "z 

- vk r ii . + 


(14) 


Here 


x 2 = 4£r + r- 


dO 2 

~cb 


1 d0 2 r 4 
r 3 dr 


(15) 


is the epicyclic frequency. For the power-law rotation Or ~ r the epicyclic frequency is simply x 2 /O 1 = 4 - q. 

In deriving these equations we have set to unity the correction factors [/?], [<3>], [Z], [£] introduced in Paper I and which take into account 
the dependence of the viscosity coefficient on temperature 77 ~ T a ' isc ( a visc = 5/2 for fully ionized gas and a visc = 1/2 for neutral gas) in 
the perturbed viscous force component N r ,N^,,N z , respectively (see Appendix A and Eq. (42) in Paper I), because the deviations of these 
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coefficients from unity have a very insignificant effect on the results. Below we shall also neglect the second derivatives with respect to z 
of the perturbed velocity components in the dynamical equations 03- 03- This can be justified if vk 2 u rj p, > v(d 2 u r ^ z /dz 2 )- To within 
a numerical factor, this inequality can be recast to the form ( k r r) 2 (zo/r ) 2 > 1. We shall see below that for k r r > 100 where the maximum 
instability increments occur and for thin discs with zo/r ~ 0.02 (see Eq. ( 149b ) this is indeed the case. 

(iii) Pressure and entropy perturbations 

In the general case by varying the equation of state Eq. © we obtain for entropy perturbations: 

Pi_ 

Po 

On the other hand, from the equation of state for ideal gas in the form p = pKT/p, we find for small temperature perturbations we have: 

El = El + Il 

Po po T 0 


N ^ Pi 

— + y— 

Cv Po 


(16) 


(17) 


(iv) Energy equation 

The linearised viscous dissipation function is 


dQ. 

Gvisc = vpr— 
dr 


dQ ~. U <P 

r - ZikrUj, - 1 — 

dr r 


+ quadratic terms. 


Here O = u^ 0 /r is the angular (Keplerian) velocity of the unperturbed flow. The linearised energy equation takes the form 

Ti 


PoKTo (. ds 0 ds 0 

- lOIJi + u z — + u r — 

p \ az or 


dQ. 

-2\k r vpnr—u ( p ■ 


■ kIczTq 


To’ 


To take into account the heat conductivity effects, it is convenient to introduce the dimensionless Prandtl number: 

vpoC p vp 0 CR/p)c p vpoCR/p) y 


Pr : 


y- 


l 


(18) 


(19) 


( 20 ) 


The Prandtl number defined by Eq. ( 120b for fully ionized hydrogen gas (y = 5/3), where the heat conduction is determined by light electrons, 
is quite low (see lSpitzeil il962h l: 


Pr, 


0.406 


1/2 


- * 0.052. 


( 21 ) 


20 • 0.4 ■ 0.225 • (2/jt)2 

Note that the presence of magnetic field in plasma decreases both electron heat conductivity and viscosity. In this case both the viscosity 
and heat conductivity are determined by ions that have larger Larmor radius than electrons, and the Prandtl number even in the case of fully 
ionized gas becomes ISpitzeil fl962h 

20 ~ 


Pi'; = 4_c p 


( 22 ) 


which is 3/8 for y = 5/3. 

In the case of cold neutral gas the Prandtl number is Pr„=2/3 accordi ng to simplifie d kinetic theory ( iHirschfelder. Curtiss & BirJ 1 9541) . 
and the heat conductivity coefficient depends on temperature as k ~ T l/2 (Spitzer 1962). 

After eliminating the temperature variations in the energy equation using Eq. ( 1161 and Eq. ( 1171 . we find 


Pi /. vk; 

- 1W + —IT 

Po \ Pr 


ds 0 ds 0 

U z ~z~ + Ur — 

oz or 


2 i k r vr(dQ./dr) p [ 

CpKTo/p U,p + po 


u) vkz 

t— + — 

y Pr 


(23) 


Here c p = ycv = y/(y ~ 1) is the specific heat capacity (per particle) at constant pressure. 

We will neglect very slow variations of the unperturbed pressure, density and entropy along the radial coordinate, i.e. set d/dr = 0 in the 
continuity equations GD, dynamic equations 03 and energy equation CD. Below we shall also denote the partial derivative with 


respect to z by prime. Thus we are left with the following system of five linearised hydrodynamic equations 
for five variable u r , u z , w^pi/po, pi/po'- 

in the anelastic approximation 

p'o 

u. — i k r u r H- u- = 0 , 

Po 

(24) 

2 • P\ 

{\o + vk/.)u r - 2Q.up = i k r — , 

Po 

(25) 

K 2 

(iw + vk 2 )u<p + —u r = 0 , 

(26) 

r , T\ Pi Po 

(\o + vkr r )u z =-H-, 

Po Po Po 

(27) 

Pi /. v ^r\ 2i k r vr(dQ/dr) 1 , 

— Ux) H-= - Uj. h- S n ll- . 

Po \ Pr / CpHTo/p c p 

(28) 
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In the right-hand side of energy equation ( 128b we have omitted the term due to pressure perturbations oc pi/p 0 because otherwise it will give 
rise to spurious unstable modes in the case of steady solid-body rotation with Q = const. 


4 DERIVATION OF THE DISPERSION EQUATION 

We start with substituting pi/po from Eq. d28b into Eq. ( 127b . Here in the right-hand side of the resulting equation two coefficients depending 
on the z-coordinate arise: 


and the Brunt-Vaisala frequency: 


©o 


Po 

— v 

Po 


■N: 


Pofo 
Po Cp 


Eil 

Po dz 




Po 


(29) 


(30) 


After differentiating the resulting equation for u z with respect to z and eliminating u'., u r and using d24b - d27b . we arrive at the following 
second-order linear differential equation for density perturbations 


Pi + Ap\ + Bpi = 0 


(31) 


with coefficients: 


A = 


HV?)' 


©o 


(dlnfl/rflnr) 


k; 


(iaj + #)(ku + vkj) - (-A?) ^ + c p (iaj + £) j 


+ — 


(itu + vk 2 ) 2 


(32) 


B = 


1 + 


1 - 


(~N 2 Z ) 


+ & n 


(cllnkl/dlnr) 


(i oj + vk 2 ) 2 
Let us introduce a new variable 


(it o+vk 2 )(icj+£) (iw + v ^) 2 Cp (icj+£) 


1 + 1 ?? 




(■ ~N ?)' 


o ’ (itu + ^)(i« + vL 2 ) - (-N 2 ) 


(33) 


P i 


Ye~^ M 


to eliminate the first derivative term in Eq. Olb : 

Y" + (B- i A 2 )Y = 0. 

This equation should be supplemented with two boundary conditions: 


Li(zo) = 0 , L,( 0 ) = 0 


(34) 


(35) 


(36) 


or 

i'zfeo) = 0, F'(0) = 0 (37) 

for two linearly independent solutions. Using the form of the coefficient A 03 it is easy to check that these boundary conditions exactly 
correspond to the physically motivated boundary conditions for pressure variations: Pi(zo) = 0 , Pi( 0 ) = 0 and pi(zo) = 0 , p\( 0 ) = 0 , 
respectively. 

This Sturm-Liouville problem 03-03 can be easily solved if coefficients A and B are independent of z. Therefore, it is necessary to 
average th ese coe fficie nts ov er som e assumed background disc vertical structure. As the background solution, we will use the polytrope discs 
with (see lKetsaris & Shakural dl998l) ): 


Po(z) = Pc 1 - — 


Po(z) = p c 


1 - \ — 

Zo 


2 \ («+l) 


T 0 (Z ) = T c 


i - I — 

zo 


(38) 


Here n is the poly trope index, p c ,p c and T c are the values of density, pressure and temperature in the disc symmetry plane, respectively. 
Therefore, the density-averaged values of quantities ©o, ©q and (—TV?) should be determined as 


£\.:)P0(Z)dZ 

where zo is the semi-thickness of the disc. Thus we obtain: 

Slar 2 

®o = (®o>-. 

Zo 


(39) 


( 40 ) 
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where, as in Paper I, we have introduced the dimensionless viscosity parameter a through the free-path length of particles Z/r and the ratio of 
the sound velocity to the unperturbed angular (Keplerian) velocity uju^ 


vk 2 

- = a(k r r) 2 , 


' U c 

a = \ — 


Uff) 


(41) 


Note that the maximum possible mean-free path in thin discs in the frame of hydrodynamic treatment should be less than the disc thickness 
Zo, i.e. Z/r = ( l/zo)(zo/r ) — (l/zo)(u s /u^) < ( uju ^). The derivative of O 0 can be written as 


The corresponding dimensionless mean values are 


®o - <®o) 


Slar 2 


<®o> = - 


n + 1 

ow<X„> 


(42) 


(43) 


<®o> = ^ (2 (n + 1)(» + 2 )B($, a visc - 1) + 


n{n + 1) 


j -2 (n+ 1 )nB(-,n- 1) 


(44) 


Similarly, for the mean Brunt-Vaisala frequency we find: 

<-^z> 

In the above formulas the dimensionless surface density of the disc is 


2 Cl 2 | 

1 n + 1 

<£oV 

1 7 


■n)B(-,n). 


r podz 

<E 0 > = - -= 2 2 "B(n + 1, n + 1), 


PcZo 


(45) 


(46) 


B(x,y) is the beta-function. When deriving these values, we have used the property that the microscopic dynamic viscosity is a function of 
temperature only voPo ~ T°' visc , with a viu: = 5/2 for fully ionized gas. For neutral gas a vjsc = 1/2 and the averaging should be performed with 
weight (see below in Section [5~2l i. 

The solution of the Sturm-Liouville problem Eq. ( 135b with boundary conditions 06} and 03 results in eigenfunctions sin(/l s z) and 
cos(/t c z) and eigenvalues T s = mn/zo and A c = (m - l/2)jt/zo (m = 1,2,...)), respectively. Substituting the eigenfunctions into Eq. ( 135b yields 
the sought for dispersion equation: 


- A 2 _ c + B - ^A 2 = 0. (47) 

The form of the coefficient A 03 immediately implies that the dispersion equation will be a tenth-order algebraic equation for a> with 
complex coefficients. It can be checked that for any value of n > 3/2 (i.e. for convectively stable disc structure), two spurious unstable modes 
arise in the case of steady-state solid-body rotation (with q = 0), which disappear if we omit the terms with (-N 2 )' (note that this term is 
not dangerous for adiabatic background structure where (-N 2 ) and its derivative vanish). Similarly unstable modes would arise if we retain 
pressure perturbations in the right-hand side of the energy equation (H3}. Due to subtleties with the energy equation in the Boussinesq and 
anelastic approximations mentioned above and analysed in JVasil et al.ll2013l) . we conclude that the Brunt-Vaisala frequency should be kept 
constant when differentiating Eq. 423. This might indicate that the background entropy gradient should be omitted in the right-hand side 
of the energy equation ( 128b from the very beginning. Not so, since its inclusion leads to the physically correct result (stabilization of the 
modes) already in the Boussinesq approximation (see Paper I). Another possibility might be the keeping both O 0 an d (~N 2 ) constant when 
deriving equation ( 13 1 } for pressure perturbations. Not so again, since the factor <T> 0 ( 129b stems from the fully leguitime linearised form of 
j-component of the Navier-Stokes equation 413. While, of course, rigorous proof of such a treatment is highly desirable, here we restrict 
ourselves to the qualitative arguments given above. Therefore, after crossing out terms with (-Air)' in Eq. ( 132b and Eq. ( 133b . we will be left 
with a sixth-order algebraic equation for a). 

The found eigenfunctions for the variable Y means that the eigenfunctions for the pressure and velocity per t urbati ons have the form 
exp(-Az/2) cos(XjZ) or exp(-Az/2) sin(Tjz), which is typical for perturbations in stratified atmospheres I Vasil et al.l 420131) . We will find that 
the maximum increment is reached for the cos mode, so substituting coefficients A and B from Eq. ( 132b and Eq. ( 133b yields the following 
quantized dispersion equation: 


(m - 1/2 )ji 

zo 


kl 


1 + - 


(ico + vk 2 ) 2 


1 - 


HY?) 


+ o' 


(<7 In kl/d In r) 


(i to + vk 2 )(iio + ) ( 1W + vk rk cJiw + ) 


-On 


(d lnO/i/lnr) 


kz 


(iaj + vk 2 ) 2 Cj>(iw+ *£) 


(i at + vk 2 ) 2 


0. 


(48) 


This is a sixth-order algebraic equation with complex coefficients (cf. cubic dispersion equation (31) from Paper I derived in the Boussi¬ 
nesq limit using local WKB-analysis). 
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5 SOLUTION OF THE DISPERSION EQUATION 


Let us analyse the solution of the dispersion equation in the anelastic approximation derived above for two fluids: the case of fully ionized 
gas with the Prandtl number Pr,, = 0.053 and Pr, = 3/8 if the magnetic field is present, and the case of neutral gas with Pr„ = 2/3. The last 
case should be treated separately, since for the dependence of the dynamical viscosity oc T 1/2 the averaging over the vertical coordinate with 
the weight po(z) is insufficient - near the surface layers of the disc the mean free-path length of particles is so large that leads to divergences 
in the term 

We will consider the laminar shear flow with the velocity profile A 2 oc r~ q so that the shear coefficient is din A/dlnr = -q/2. All 
relevant frequencies will be normalized to the local Keplerian value A and denoted with tilde. In the numerics the adiabatic index of the gas 
is set to y = 5/3. 

It is convenient to write down the dimensionless dispersion equation for the dimensionless mode frequency ai as a function of the 
dimensionless variable ( k r r ) with dimensionless viscosity parameter a and the disc thickness 


£o 

r 



(49) 


Here the dimensionless factor Hi takes into account the vertical disc structure, and in the case of the polytrope accretion discs ^ = 2{n + 1) 
iKetsaris & Shakurail 19981) . The values of the sound velocity u s then should be taken in the disc symmetry plane. 

The dispersion equation (48} in the dimensionless form reads: 


(m - l/2)it\ 2 

zo/r / 


(k r r) 2 


1 + 


(i< 2 > + a(k r r) 2 ) 2 


<®o )a 


(-N 2 ) 


<®o )a 


(-q/2) 


(id) + a(k r r) 2 )( id) + ° ( p r) ) (zo/ r ) 2 (i® + a(k r r)~Y Cp (ia> + ) 


( ~q/2)a 


tK& 


(zo/r) (id) + a(k r r) 2 ) 2 Cp(i£) + 2) / + _£ 


(id) + a(k r r) 2 ) 2 


= 0 . 


(50) 


5.1 Fully ionized gas 

The results of the solution of the dispersion equation El are shown in Fig. |T| for two background vertical structures of a thin Keplerian 
disc. In the left panel of Fig.Qjwe present the solution for the polytrope disc structure with constant entropy described by the polytrope index 
n = |. Here the Brunt-Vasala frequency N 2 vanishes. It is seen that the unstable anelastic mode (the one with the negative imaginary part in 
the bottom panel) arises at k r ~ 50 - 150, where the approximation of the incompressibility ( dp/dt = 0) is justified. In fact, this is two modes 
with equal absolute values but different by sign real parts that demonstrate the overstability (marked with arrows in the upper left panel). This 
is different from the Boussinesq limit considered in Paper I, where one of the Rayleigh inertial modes became viscously overstable in the 
presence of viscosity. In the anelastic approximation, the Rayleigh inertial modes remain always stable, and the unstable modes are split out 
from the Rayleigh modes. In the right panels of Fig.[T]we present the solutions for the background disc structure with vertically increasing 
entropy (shown is the solution for n = 2), which is convectively stable (N: > 0) in the absence of viscosity and shear. However, the presence 
of even small viscosity makes the Keplerian flow convectively unstable even in this case. 

In Fig. [ 2 ] we explore the effect of different dimensionless parameters of the problem on the increment of the overstability. First, in the 
left panel of Fig. [2] we study the effect of changing the Prandtl number, which describes the dumping effect of thermal conductivity on small 
perturbations. For fully ionized gas, the Prandtl number is maximal when the small (but still dynamically unimportant) magnetic field is 
present, and both the viscosity and heat conductivity are mediated by ions which have larger Larmor radius than electrons (Pr,=3/8, see 
above). We also show the results of the dumping effect of possible radiative conductivity parameterized in terms of the effective Prandtl 
number (Pr/2 and Pr/11, see Eq. (47) in Paper I). The smaller the Prandtl number, the smaller the instability increment, which is physically 
clear. Second, in the central panel of Fig.|2]we show the effect of changing the viscosity (parametrized in terms of the effective mean free-path 
length of particles, l/r). The larger the viscosity, the higher the instability increment. Finally, in the left panel of Fig.[2]we demonstrate the 
effect of the disc semi-thickness (parametrized by the central sound speed to the angular velocity ratio, uj up. At a given mean-free path of 
particles, the thinner the disc, the higher the increment. 


5.2 Neutral gas 

For fully neutral gas with the Prandtl number Pr„ = 2/3 we should first make new averaging of quantities ®o and (—/V?) with weight Pq(z) 

to avoid divergence of the value Of, near the disc surface due to very large mean-free path length of particles in the polytropic discs. We find: 


«®o» = - 


n + 1 

(&visc "t" U)(/\Yp ) 


(51) 


«®o» = 


«£o» 


2 (n + 1 )(n + 2 )B(-,a vist 


1 + n/2) + 


n{n +1) 

Q-visc 1 V II 


■ 2(n + 1 )nB(—,n - 1 + n/2) 


(52) 
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Figure 1. Left: Real and imaginary parts of anelastic modes in fully ionized gas with electron heat conductivity (Pr^ =0.052) and viscosity parameters u s ju,p = 
0.01, Z/r = 10 -4 for the adiabatic density distribution (n = |). Right: The same for vertically increasing entropy distribution with n = 2. 


q=3, n=3/2, u/u^O.OI, l/r=1 (T 4 q=3, n=2, uAj^O.OI, Pr e =0.052 q=3 , n=2j ,/ r=10 -4 r^o.052 





Figure 2. Left: Imaginary part of the viscously unstable anelastic mode in fully ionized gas with different Prandtl numbers. Viscosity parameters Ug/u^ = 
0.01, l/r = 10 -4 for vertically increasing entropy distribution with n = 2. Middle: The same for different mean free-path length of particles l/r with fixed disc 
thickness parameter u s /U(p = 0.01. Right: The same with different disc thickness parameter u s /u,p for fixed mean free-path length of particles l/r = 10 -4 . 
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Figure 3. Left: Real and imaginary parts of anelastic modes in neutral gas for the adiabatic entropy distribution (n = |). Right: The same for vertically 
increasing entropy distribution with n = 3. Note that the real part of the solution is quite insensitive to the poly tropic structure of the disc. 


, 2 ln+ 1 \ 3 

<<( -^ )>>= (^>>(—-T ( 5’ 2 ' 7) - (53) 

«X 0 » = 2 4, 'S(2« + 1,2« + 1), (54) 

In Fig. [3]we show the real and imaginary parts of anelastic modes for the case of neutral gas with the Prandtl number Pr„ = 2/3. The 
disc thickness is zo/r ~ = 0.01, the mean free path of ions is l/r = 1CF 4 . Since the Prandtl number is quite large, unlike the case of 

the ionized gas shown in Fig.|T| the instability is present even for the quite significant vertical background entropy gradient ( n = 3, the right 
panel of this Figure). 


6 DISCUSSION AND CONCLUSIONS 

In the present paper we have extended the modal analysis of small axially symmetric perturbation in sheared Keplerian flows with micro¬ 
physical viscosity and heat conductivity, which we started in Paper I. In contrast to Paper I, where the Boussinesq approximation was used, 
here we have formulated the problem in the anelastic approximation and taken into account vertical boundary conditions in thin discs. In this 
approximation we have obtained the second-order linear differential equation with respect to the z-coordinate for small pressure perturba¬ 
tions, Eq. Oil , with coefficients (Eq. ( 132b and Eq. ( 1331) ) depending on the height above the disc plane z. We have assumed the background 
polytropic vertical disc structure and averaged the coefficients Eq. ( 132b and Eq. ( 1331 with weight p 0 (z) over the vertical disc height. This 
allowed us to solve the Sturm-Liouville problem to obtain the discrete spectrum of eigenvalues and eigenfunctions (cos(T c z) or sin(/ljz)). 
Substitution of these functions into Eq. ED resulted in the dispersion equation for normal modes of axially symmetric perturbations along 
the radial coordinate, Eq. (48]. This turned out to be an algebraic sixth-order equation, solution of which for a wide range of the viscosity 
parameters in thin Keplerian discs are presented in Fig. l 1 131 Note that in the Boussinesq limit the dispersion equation had only the third order. 
We have found that in a wide range of wavenumbers k r r » 1 two unstable modes are split from the classical inertial Rayleigh modes (in the 
Boussinesq limit one of the Rayleigh modes displayed the overstability in the presence of viscosity). Qualitatively, the results of the present 
paper are in agreement with findings of Paper I. However, in contrast to Paper I, where the Boussinesq approximation was used, the results 
are quantitatively different. 
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The description of the hydrodynamic flows in the anelastic approximation, although neglects the sound modes, is more precise and 
takes into account the important term in the continuity equation, (l/p 0 dpo/dz)u-. This means that the analysed small perturbations are no 
more purely transversal, as in the Boussinesq limit. However, it is important to note that in both anelastic and Boussinesq approximations 
the pressure variations pi/po should be neglected in the energy equation, otherwise fictitious unstable solutions emerge in the case of the 
steady-state solid-body rotation. We have also found that the overstability appears in the cases where there is a non-zero vertical gradient of 
the quantity (p' 0 /p 0 )v ~ -(Q. 2 z/T 0 )v. 

In both approximation we have found that the increment of overstable modes increases with viscosity and the background vertical 
pressure gradient, suggesting a convective nature of the overstability: the viscous heat generation in the sheared flow in the gravitational field 
around a central star makes the flow convectively unstable. It is tempting to suggest that this instability may be the seed for turbulence in 
Keplerian discs even in the absence of magnetic fields. 

Note the many faces of the viscosity in the considered problem. The higher the viscosity in the right-hand side of energy equation 
(128}. the stronger the viscous energy generation due to the shear leading to the buoyancy of the perturbed regions. On the other hand, in 
the dynamic equations {ED- G3 the viscosity damps the perturbations. In the right-hand side of these equations, we have neglected terms 
vu ", vu'^,vu" with second derivatives of perturbed velocities. We expect that their taking into account will somewhat decrease the increment 
of the viscous-convective instability and narrow the interval of unstable wavenumbers in Figs.|T]and[3] If we retain these second derivatives, 
we will obtain a much more complicated sixth-order differential equation for perturbations. This is a separate problem to be addressed 
elsewhere. Here we have restricted ourselves to solving only the second-order differential equation Eq. d3 1 b . 
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